#ifndef __line_hh__
#define __line_hh__

#include "nodal_fe.hh"
#include "fem_traits.hh"

#include <iostream>

using namespace std;

/// Class for constant FE on interval
template <unsigned int ORDER, fastfem::PolynomialFamily PFAMILY, typename T=double>
class SegmentFiniteElement : public NodalFiniteElement<1, fastfem::Segment, ORDER, PFAMILY, T> 
{
public:
   SegmentFiniteElement() : NodalFiniteElement<1, fastfem::Segment, ORDER, PFAMILY, T> () {}
   virtual ~SegmentFiniteElement() {}

  /** virtual function which evaluates the values of all
       shape functions at a given point ip */
  virtual T evaluate_basis ( unsigned int , const Point<SPATIAL_DIMENSION, T>& )
  {
    femlib_assert(false, "Not implemented");    // throw an error deliberately
    return T();
  }
  
  /** virtual function which evaluates the values of all
       shape functions at all point ips */
  virtual void evaluate_basis_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<T>& basis_vals )
  {
    femlib_assert(false, "Not implemented");    // throw an error deliberately
  }
  
  virtual void evaluate_gradient ( unsigned int ib, const Point<SPATIAL_DIMENSION, T>& pt, Point<SPATIAL_DIMENSION, T>& grad)
  {
    femlib_assert(false, "Not implemented");    // throw an error deliberately
  }
  
  virtual void evaluate_gradient_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<Point<SPATIAL_DIMENSION, T> >& grad_vals )
  {
    femlib_assert(false, "Not implemented");    // throw an error deliberately
  }
  
  virtual void evaluate_hessian ( unsigned int ib, const Point<SPATIAL_DIMENSION, T>& pt, PointList<SPATIAL_DIMENSION, T>& hess)
  {
    femlib_assert(false, "Not implemented");    // throw an error deliberately
  }
  
  virtual void evaluate_hessian_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<PointList<SPATIAL_DIMENSION, T> >& hess_vals )
  {
    femlib_assert(false, "Not implemented");    // throw an error deliberately
  }
  
  void reference_element_coords(std::vector<Point<SPATIAL_DIMENSION, T> >& nodes)
  {
    femlib_assert(false, "Not implemented");    // throw an error deliberately
  }
  
  virtual int test()
  {
    std::vector<Point<SPATIAL_DIMENSION, T> > nodes ;
    this->reference_element_coords(nodes);
    
    cout << "\n -- SegmentFiniteElement --";
    cout << "\n   Number of reference element nodes:= " << nodes.size() ;
    for (unsigned int i=0; i < nodes.size(); ++i)
      cout << "\n        nodes[" << i << "] = [" << nodes[i][0]  << ", " << nodes[i][1]  << ", " << nodes[i][2] << "]" ;
    cout << endl ;
    return 0 ;
  }
  
};


template <unsigned int ORDER, fastfem::PolynomialFamily PFAMILY, typename T>
struct FEMTraits<SegmentFiniteElement<ORDER, PFAMILY, T> >
{
  
  enum { FE_Dimension = 1 } ;
  enum { FE_Family = PFAMILY } ;
  enum { FE_Order = ORDER } ;
  enum { N_Nodes = ORDER+1 } ;

  typedef T BaseType;
  typedef SegmentFiniteElement<ORDER, PFAMILY, T> FEType;
  typedef NodalFiniteElement<1, fastfem::Segment, ORDER, PFAMILY> ParentFEType;

};


// ALL Lagrange type 1-D basis functions: Linear elements

/** virtual function which evaluates the values of all
      shape functions at a given point ip */
template <>
inline
double SegmentFiniteElement<1, fastfem::Lagrange, double>::evaluate_basis (  unsigned int ib, const Point<SPATIAL_DIMENSION>& px)
{
#define NPTS 2
  femlib_assert(ib < NPTS, "Invalid basis index");

  switch(ib)
  {
    case 0:
      return (1.0 - px.x()) * 0.5 ;
    case 1:
      return (1.0 + px.x()) * 0.5 ;
    default:
      cout << "\nInvalid basis index. Only 2 available.\n";
      exit(2);
  }
#undef NPTS
}

/** virtual function which evaluates the values of all
      shape functions at all point ips */
template<>
inline
void SegmentFiniteElement<1, fastfem::Lagrange, double>::evaluate_basis_all ( const PointList<SPATIAL_DIMENSION>& ips, std::vector<double>& basis_vals )
{
#define NPTS 2
  basis_vals.resize(NPTS*ips.n_points);
  for (unsigned int i=0, j=0; i < ips.n_points; ++i, j+=NPTS) 
  {
    basis_vals[j] = (1.0 - ips.x(i)) * 0.5 ;
    basis_vals[j+1] = (1.0 + ips.x(i)) * 0.5 ;
  }
#undef NPTS

  return ;
}

template <>
inline
void SegmentFiniteElement<1, fastfem::Lagrange, double>::evaluate_gradient (  unsigned int ib, const Point<SPATIAL_DIMENSION>&, Point<SPATIAL_DIMENSION>& grad)
{
#define NPTS 2
  femlib_assert(ib < NPTS, "Invalid basis index");

  switch(ib)
  {
    case 0:
      grad[0] = -0.5 ;
      break;
    case 1:
      grad[0] = 0.5 ;
      break;
    default:
      cout << "\nInvalid basis index. Only 2 available.\n";
      exit(2);
  }
#undef NPTS
}

template <>
inline
void SegmentFiniteElement<1, fastfem::Lagrange, double>::evaluate_gradient_all ( const PointList<SPATIAL_DIMENSION>& ips, std::vector<Point<SPATIAL_DIMENSION> >& grad_vals )
{
#define NPTS 2
  grad_vals.resize(NPTS*ips.n_points);
  for (unsigned int i=0, j=0; i < ips.n_points; ++i, j+=NPTS) 
  {
    grad_vals[j][0] = -0.5 ;
    grad_vals[j+1][0] = 0.5 ;
  }
#undef NPTS

  return ;
}

template <>
inline
void SegmentFiniteElement<1, fastfem::Lagrange, double>::reference_element_coords(std::vector<Point<SPATIAL_DIMENSION> >& nodes)
{
#define NPTS 2
  nodes.resize(NPTS);
  nodes[0][0] = -1.0 ;
  nodes[1][0] = 1.0 ;
#undef NPTS
}


// ALL Lagrange type 1-D basis functions: Quadratic elements

/** virtual function which evaluates the values of all
      shape functions at a given point ip */
template <>
inline
double SegmentFiniteElement<2, fastfem::Lagrange, double>::evaluate_basis (  unsigned int ib, const Point<SPATIAL_DIMENSION>& px)
{
#define NPTS 3
  femlib_assert(ib < NPTS, "Invalid basis index");

  double x = px.x();

  switch(ib)
  {
    case 0:
      return x * (x - 1) * 0.5 ;
    case 1:
      return x * (1.0 + x) * 0.5 ;
    case 2:
      return (1.0 - x * x) ;
    default:
      cout << "\nInvalid basis index. Only 3 available.\n";
      exit(2);
  }
#undef NPTS
}

/** virtual function which evaluates the values of all
      shape functions at all point ips */
template<>
inline
void SegmentFiniteElement<2, fastfem::Lagrange, double>::evaluate_basis_all ( const PointList<SPATIAL_DIMENSION>& ips, std::vector<double>& basis_vals )
{
#define NPTS 3
  basis_vals.resize(NPTS*ips.n_points);
  for (unsigned int i=0, j=0; i < ips.n_points; ++i, j+=NPTS) 
  {
    const double x = ips.x(i);

    basis_vals[j] = x * (x - 1) * 0.5 ;
    basis_vals[j+1] = x * (1.0 + x) * 0.5 ;
    basis_vals[j+2] = (1.0 - x * x) ;
  }
#undef NPTS
  return ;
}

template <>
inline
void SegmentFiniteElement<2, fastfem::Lagrange, double>::evaluate_gradient (  unsigned int ib, const Point<SPATIAL_DIMENSION>& px, Point<SPATIAL_DIMENSION>& grad)
{
#define NPTS 3
  femlib_assert(ib < NPTS, "Invalid basis index");

  switch(ib)
  {
    case 0:
      grad[0] = px.x() - 0.5 ;
      break;
    case 1:
      grad[0] = px.x() + 0.5 ;
      break;
    case 2:
      grad[0] = - 2.0 * px.x() ;
      break;
    default:
      cout << "\nInvalid basis index. Only 3 available.\n";
      exit(2);
  }
#undef NPTS
}

template <>
inline
void SegmentFiniteElement<2, fastfem::Lagrange, double>::evaluate_gradient_all ( const PointList<SPATIAL_DIMENSION>& ips, std::vector<Point<SPATIAL_DIMENSION> >& grad_vals )
{
#define NPTS 3
  grad_vals.resize(NPTS*ips.n_points);
  for (unsigned int i=0, j=0; i < ips.n_points; ++i, j+=NPTS) 
  {
    const double x = ips.x(i);

    grad_vals[j][0] = x - 0.5 ;
    grad_vals[j+1][0] = x + 0.5 ;
    grad_vals[j+2][0] = - 2.0 * x ;
  }

#undef NPTS
  return ;
}

template <>
inline
void SegmentFiniteElement<2, fastfem::Lagrange, double>::reference_element_coords(std::vector<Point<SPATIAL_DIMENSION> >& nodes)
{
  nodes.resize(3);
  nodes[0][0] = -1.0 ;
  nodes[1][0] = 1.0 ;
  nodes[2][0] = 0.0 ;
}



/*

inline
template <unsigned int order>
void SegmentFiniteElement::compute_volume()
{
  volume = fabs(vertices[3] - vertices[0]) ;
  
  // check if vertices are degenerate
  if ( volume < 1e-8 )
  {
    cout << "\n";
    cout << "Line - Fatal error!\n";
    cout << "  Element has zero volume.\n";
    cout << "  Area = " << volume << "\n";
    cout << "\n";
    exit ( 1 );
  }
}

inline
double Line::element_volume()
{
  return this->volume ;
}
*/

#endif // #ifndef __line_hh__
